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decomposition method. PADOC comes along with a new solution approach, namely, 
average NonLinear Programming (aNLP). The aNLP theorem suggests generating 
staircase optimal solution (low resolution) and turning it into a distributed optimal 
solution (high resolution) by exploiting the Hamiltonian. The analytic architecture 
promotes PADOC to institute a certain analytic connection between the discrete 
nodes. This renders stability, accuracy within an analytic resolution (for the mathe- 
matical objects), robustness and a cut-down in the number of the decision variables 
in the frame of NLP. We also give a rigorous proof that the multipliers associated 
with Karush-Kuhn- Tucker optimality conditions in the frame of NLP, as transcribed 
by PADOC, correspond to a backward analytic solution for the costate equations. We 
finally let PADOC speak through some examples. 
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1 | INTRODUCTION 


The XXI century has witnessed innumerable applications of optimal control theory in disparate areas of research such as 
biomedicine, e. glil aircrafts/spacecrafts, e. g 201801 and robotics, e. gE In this respect, the associated industries are in demand 
of solution methodologies for Optimal Control Problems (OCP) with high accuracy, robustness and reasonable computational 
time. In a conventional point of view, the methods for OCPs can be put into two categories, indirect and direct approaches"! 
The indirect approach is linked to Pontryagin Maximum Principle (PMP) leading to a system of Multi-Point Boundary Value 
Problem (MPBVP). In a direct approach, the transcribed/parametric version of the OCP as a NonLinear Programming (NLP) is 
fed to an optimization solver such as interior point or Successive Quadrative Programming (SQP). Indeed, the method of tran- 
scription affects the sparsity of the Karush-Kuhn- Tucker (KKT) system and accordingly the search direction matrix, accuracy 
and the required iterations in the frame of NLP. However, as long as the transcription method meets a sufficient level of accuracy 
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over the search space, the general structure of the NLP problem is preserved, that is, a convex OCP remains convex on dis- 
similar accurate transcription methods. In this respect, a rigorous mathematical analysis is still lacking to evaluate the possible 
repercussions on the NLP due to different transcription methods (a more detailed review on this matter can be tracked inb). 

The advantage of the indirect approach is to provide accurate solutions. However, theoretical analysis reveals a mathematical 
connection/equivalence between the indirect and direct methods?l Generally speaking, Covector Mapping Principle (CMP) 
asserts that there exists a linkage between an indirect and a direct method through Hahn-Banach theorem, see e. g. Lola andli3l 
In other words, CMP (which is an order-preserving map between the primal and dual variables) pairs a direct method with an 
equivalent indirect method. This means, the solution by a direct method converges to the solution by its indirect counterpart. 
This understanding has led to a newly coined term, Hamiltonian programming (for more information, interested readers are 
referred tol and the relevant references therein). Moreover, CMP is regarded as a criterion to evaluate whether a transcription 
method is complete?l 

The present research aims to propose a novel analytic transcription approach for optimal control problems that is complete 
(accounting CMP criterion). The corestone of the presented method, Perturbed Analytic Direct transcription for Optimal Control 
(PADOC), is to provide a rigorous analytic connection between the intermediate nodes. In this sense, solution for the state 
variables is exact. The analytic connection can be explained by assuming two adjacent nodes from the state dynamics with one 
decision variable being allowed to vary in parameter optimization. In this situation, PADOC assures that the two nodes are 
always analytically connected over the entire search space. Therefore, PADOC is guaranteed for the convergence of the system 
dynamics with only a few segments for a wide class of nonlinearity (the property lacking in most of the numerical integration 
methods). Comparing to the collocation methods, the immediate understanding is to provide a transformed cost functional within 
an analytic resolution. In the sequel, we provide a state of the art critique on transcription approaches to highlight the existing 
gaps. 

A direct trajectory optimization requires a transcription method to transform the original problem into an NLP. The tran- 
scription methods can be viewed as shooting and collocation. А circumstantial review of these methods can be found in several 
resources (see e.g. ПЭПб) and this suffices to avoid superfluous description. Specifically, shooting techniques employ either 
explicit (such as Euler and Runge- Kutta) or implicit (such as Hermite-Simpson) approaches to discretize the system dynamics 
normally with piecewise constant control (see e. g 7 and 8h, In collocation techniques (see e. в.П9100 and, the system dynam- 
ics and control are considered as piecewise functions, approximated by orthogonal polynomials, and the system dynamics are 
satisfied at some nodes, known as the collocation points. This leads to an NLP with states being explicitly dependent on the 
control resulting in more memory storage’ As long as the stability of the solution to the system dynamics is concerned, col- 
location methods are more advantageous over the shooting methods. This statement is supported by referring to the fact that 
in collocation methods, the dynamics’ path between the imposed handling nodes are governed/interpolated by some predefined 
base functions of polynomial types. This is while shooting methods need to estimate the dynamics’ path by recursive numeri- 
cal drifts. Therefore, stability of the shooting methods is restricted to the stability of the numerical discretization. On the other 
hand, solutions by the collocation methods are valid only at the collocation points. A drawback of the shooting methods is the 
need for a prior knowledge on the sufficient discrete segments to guarantee the convergence of the system dynamics. Therefore, 
the number of grids is normally considered large enough to guarantee convergence, leading to extra memory storage. Besides, 
the sole concept of convergence, apart from the accuracy, can be loosely interpreted in such a way that the convexity of the 
OCP is preserved. However, the Hessian search direction matrix is altered affecting the required number of Newton-type iter- 
ations to find the optimal solution with a limited accuracy. One of the underlying reasons associated with the failure of the 
optimization module is that a numerical discretization cannot guarantee that the system dynamics are satisfied in iterative steps 
performed by the optimization module. In the case of singular OCPs, the Hessian matrices are ill-posed. Experience shows that 
this ill-conditioning is magnified as the number of nodes increases. This is whilst, there are many practical nonlinear OCPs 
requiring a large number of nodes in order to secure the convergence of the system dynamics mostly in the frame of the numer- 
ical shooting transcriptions. An amalgamation of these statements indicates a plausible failure of the shooting methods in these 
circumstances. Therefore, both collocation and shooting methods present well-identified drawbacks that we claim to (at least 
partially) overcome with PADOC for different classes of OCPs. 

A more challenging situation is the presence of singular arcs in the OCP, leading to the ill-posedness of the Hessian matrix. In 
other words, the projected Hessian matrix is not positive definite22l This ill-conditioning shows itself as an oscillatory or fluc- 
tuating response in the singular region. There are some reported strategies to overcome this situation, e.g., mesh refinement and 
hubristic strategies (see e. g. ZLAB, combining indirect and direct pathways (see e. в.29129180)), parameter continuation 
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(sometimes referred as homotopy) or a combination of these strategies 2l More specifically 3l suggests a mesh refinement strat- 
egy to increase the resolution of a singular arc solution in a direct method. The nested direct transcription optimization as in26] 
splits the OCP into an inner and an outer OCP and finds an optimized mesh distribution. 1220] the authors proposed a strategy 
that amalgamates indirect and direct methods. More specifically, the OCP is initially solved by direct optimization to extract 
the solution structure and to estimate the onsets of the costate variables. When it comes to parameter continuation? the strat- 
egy is to switch a singular OCP to a nonsingular OCP. Basically, parameter continuation consists of smooth deformations from 
a known (or easy to solve) OCP to the problem of interest in infinitesimal steps such that the solution of each step serves as 
a meaningful initial guess for the next step. In this way, parameter continuation is recommended as a refinement strategy in a 
direct or an indirect method.22l contains mesh refinement processing as well as a regularization procedure with an embedding 
parameter to avoid ill-conditioning of the Hessian matrix by virtually turning a singular problem into a nonsingular one. 

In a recent publication? an analytic transcription method presumably as the first appearing attempt in the literature for 
trajectory optimization problems has been presented. The method has been launched upon Homotopy Analysis Method (HAM) 
that is a powerful analytic nonlinear solver over about two decades (interested readers are referred to some original references 
such as (see e. g. BALBB), However, it certainly cannot guarantee a global convergence for the nonlinearities. It is due to the 
fact that HAM normally seeks a unified solution for the entire domain of interest. Therefore, the convergence of HAM is reliant 
on the existence of a fixed point for the inverse of the linear operator as it applies to the nonlinearity over the entire domain. 
Since the choice for the linear operator is free in HAM and there is no theory to restrict this choice, it does not guarantee a 
global convergence. In this respect, authors have treated different nonlinearities with different linear operators through HAM 
on a case by case basis (see e.g. B466] and33), Therefore, the method proposed іп will fail as it confronts strongly-nonlinear 
system dynamics. As a consequence, PADOC stands as the first robust analytic transcription method in the literature. 

We resort now our discussion to the assessment of the completeness of PADOC using the CMP principle?! Connection 
between the multipliers associated with Karush-Kuhn- Tucker (KKT) optimality conditions and the continuous costate functions 
has been investigated in the past in several resources such as39]H0l angl] For some transcription methods, such as the ones 139] 
апа the multipliers are connected to a discrete numerical solution for the costate variables (dual variables) with a reduced 
order of accuracy compared to the accuracy of the state variables (primal variables). For the Legendre-Gauss-Lobatto (LGL) 
collocation, as in“! the costates at the LGL points are equal to the KKT multipliers divided by the LGL weights. Here, we prove 
that the KKT multipliers at the SADM nodes correspond to a backward analytic solution of the costate equations. Moreover, we 
prove that the cost functional in the frame of NLP is decomposed in a pure analytic manner. 

Therefore, the novel contributions of the present paper can be viewed as: 1) presenting SADM as a new analytic nonlinear 
solver and 2) a new NLP solving strategy, namely average NonLinear Programming (aNLP). SADM and aNLP are exploited to 
form a novel analytic transcription method, PADOC. 

PADOC has the following attributes: a) to be hp adaptive. b) to have the minimum number of nodes. c) to show marginal 
sensitivity to the initial guess due to the analytic connection between the imposed nodes. d) to show highly accurate estimations 
for the three main mathematical objects, namely, the system dynamics, constraints and the cost functional. e) to be capable 
of handling singular OCPs. f) to exhibit a quite spare programming for regular OCPs (no bang/singular structures). g) to have 
primal-dual consistency (CMP completeness). 

Finally, itis worth mentioning that, in general, PADOC could be launched upon other analytic nonlinear solvers such as HAM. 
However, for the sake of convergence, it is a mandate to design a segmented form of HAM. Even though the usage of HAM 
within PADOC lays outside the present paper. We will aim in the future to present different forms of PADOC with possible 
distinguished merits. 

The paper is constructed as follows: In SecD]we present ADM, as a prominent analytic nonlinear solver. In SecB] we design 
SADM as anew analytic solver together with the supportive theorems[I]and]to assert the unconditionally-convergent argument. 
Sec/4] is devoted to the mathematical representation of SADM in OCPs. In Sec[5] we present the primal-dual consistency of 
PADOC together with the aNLP theorem. In Sec. [6] the PADOC algorithm is given. In Sec[7] we solve some examples with 
varying complexity to show the applicability and robustness of PADOC. These examples are: 1) an orbital transfer problem 
(7.1); 2) ће Van der Pol oscillator (7.2); 3) Brachistochrone problem (73); 4) Aly problem (7.4); 5) maximum range of a hang 
glider (7.5); 6) minimum time aircraft trajectory in climbing phase (7.6). 
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2 | ADOMIAN DECOMPOSITION METHOD 


Let us consider a nonlinear problem in a general form as , 
LO] = V [x] +®[х@] +8) teo cR (1) 


k 
In above, £ is an invertible linear operator, e.g., £ = =( . ) with k being the highest order of derivative, R is the remainder 
of the linear operator, M is the nonlinear operator, and g(t) is a known source function. We point out that the choice for the 
linear operator is not unique. For example, with regard to Lane-Emden equation, (х@)) + 2 + £ (xq) + g(x) = 0, one сап 
2 
choose any of the following linear operators: £ = =( . ) Or 04 (12 (9). 
Eq. (1) is accompanied by a set of boundary conditions as , 


xO) = у(00, 270) =70Q), xO) = 72), n О лу yk - 1) (2) 


In order to restrict the statement of the problem to ГУР] it is assumed that , 


lysi ТЕЖЕ Е (3) 


In principle, ADM considers an inverse operator as, 


mere (4) 


The inverse operator transforms Eq.(1) into , 
k-1 


d 1 d? 2 l (p 
x(t) — x(t) — 27000 — to) = БП ag ot =) = a Pig) — ta)? (5) 
= L! (N [x] RED) + £L (g0) 
The solution and the nonlinearity are decomposed as, 
х@ = AHOA — Mix, N [xO] = }ў,.А,(х@)) (6) 
i=0 i=0 
In above, А, are Adomian terms. To obtain the Adomian terms, we expand N [xc] about the function x1), 
d (х= xo)? d? (х= хо) 43 
N [xc] = N [xo] + (х – xo^ bo] + ar ga bed + Еа [xo] Ais (7) 
Therefore, 
со 2 со 
V d (еа л) а? LAE. i) d? 
N [xc] = N [xo] + (2) aN bel + — ам! o] + 3! ax" Pe ee (8) 
After some manipulations, we obtain , 
со k-n 
1 
N [xc] = N [xo] + by (2; di > Ön cp and Xx JI [xo] ) (9) 
n=1 k=l “`7 б, 
Where,ó;. jis Kronecker delta , 
ZEN ао) 
| 0 izj 


It is worth mentioning that ADM can also be adjusted to Boundary Value Problems (BVP). This adjustment is known as Duan-Rach (D-R) ADM (secl апа). То 
track a frontier nonlinear decomposition solver, we refer to Optimal Decomposition Method (ODM) recently proposed by Jafarimoghaddam, Soler et a1 44 
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Eq. (9) together with Eq.(6) gives , 


k=n 
1 
А, = x xi X б, i kis. А, gu) [xo] 
kad dubi, 
Therefore, 
Ao(xq(t)) = N (xo) 
у 
2 
а хү d? 
As (xo), x1, 00) = ху n) + ara o 
d d xi 
Aa (xo(D. x1 (0, х200), хз(0)) = эл) + к 09) а N (xo) 
Adomian and Rach proposed® 
1 d” so " 
А, (хо, xi (D. х0) = = 2 [NV ( 2 хур!)|, 


Connecting Eq. (5) and Eq. (6). the recursion scheme becomes, 


ld 


k-1 
X(t) = x(tg) + < 000 — to) + 2192 x(tg)(t — t9) + ЮУ хх” — to) + £^ (g(t) 


х„@ = L (R[x,1]) - £ ! (A4), 221 


(1) 


(12) 


(13) 


(14) 


ADM can be shown to diverge in many practical nonlinearities (see e.g., SecB.1). To circumvent it, we propose a modified 


version (SADM). 


3 | SEGMENTED ADOMIAN DECOMPOSITION METHOD (SADM) 
In SADM, the domain of interest is divided into some intervals, 
N 
Q= U Q;(t)-14)) 
j=l 


The recursion scheme in the first interval is obtained by, 


k-1 
o, O = хб) + D Fx Pr to)” + La! (0) 
р=1 `` 


Xo O= £a (R[x,-1,0,1) + £a (А, 1), п> 1 


Therefore, the solution in the first interval is, 
n 
Xo, (0 = Хо,о, + X19, Fas = * Хо, 
= 
For the proceeding intervals, 
k-1 1 
= (p) р -1 
хоо (0) = xo, (5) + У тха, 07-004 =t) + £g (s) 
p=1 *` 


х„о,@ = £g (Ж[х„о]) + £g (Ani), п>1 
Thus, 


n 
Xo (f) = хоо, + хо, +... = È хо,» jz2 
= 


(15) 


(16) 


(17) 


(18) 


(19) 
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The Adomian terms are computed by, 


1 d" = | 
А, (хоо, (0), хо (D. x49 (0) = a T. [N( > Xp), J21 Q0) 
d i=0 


x(tg) and x? (tg), р = 1,2, ..., are known form the boundary conditions and n is the order of perturbed fields in each interval. 
Theorem 1. The analytic series solution by SADM is convergent on Q, if Bp is a contraction , that is, 1517 || < 1. 


Proof: We can reconstruct the nonlinearity over the interval Q, as, 


£g [xo] - No, [x0] =e, NM=N+R (21) 
Since, M is a nonlinear Operator, so is N . Therefore, 
(Ухо) - DA =в@) (22) 
i=0 i=0 
From Eq. (22). 
Хоо = Q(t) 
X29; = аА! (23) 
2p 
Хо, = Q; Aa 


In above, from Eq. (18). Gt) = Bt) + b g(t), where, /3(7) is obtained from the boundary conditions. 
For a truncated series, we can write , f 


NE= Y digas). =) tas S oa Q4) 
i=0 i=0 i=l 

Therefore, 

хто = $— £g А (хоо) = СА (xoa, ) 
хо, + Хо, = 5 = Ly Ay (хоо) + £s A, (хоо X19,) = £5 A (хо, t xio) 
(25) 
хо, +%2,0, ++ Хо, = Sua = СА (хоо) + £g A, (xoa, Xia) + + 
£o Â, (Xon, X10, +» Xna) = £o X (хоо, *5,) 

The recurrent relation, 5,, = p (хор, + S, ),сап be written as, S = EN (хов, +S ). Therefore, S is the fixed point 
of £N and the sequence is conelgent i.e., lim, ,,, (S — 5,) > Oif ca is a contraction , that is, [ I < 1. 
Theorem 2. Let II = ү[°,] = X I i]. ] € j € N bethe Banach space. Assume that 7 (t, x) = £3 N [x] satisfies Lipschitz 
condition that is, ||7 (t, u) — T(t, v)|| € x|ļu —v|]j Vu, v € II, then, Eq. has only one solution if к<1 
Proof: let и and v be two distinct solutions of Eq. Q1). Then we can write, 

и= T (tu) + B, + Lg 8) 


(26) 
= T (t, v) + 8,0) + £g g() 
since /3, (1) and B,(t) are obtained from the boundary conditions, we have, З (т) = B,(t). Therefore, 
llu — oll = 1170, u) — TG, 0) < lu — vll Q7) 


Since, к < 1, then, u = о. 
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Corollary 1. The operator LN is not necessarily a contraction. This restrains the applications of the standard ADM. In 
other words, the standard ADM is convergent if a fixed point theorem holds for LN over the entire domain. A fixed point 
theorem exists for continuous monotonic functions such as a corollary on Knaster-Tarski fixed point theorem allocating it to the 
case of monotonic functions and Kleene’s theorem with the functions being additionally continuous (secl and), However, 
the existence of a fixed point theorem does not hold for a general non-monotonic function except for some specific classes. 
Therefore, if the solution to the nonlinearity can be assumed as a combination of segmented monotonic functions, such that 
15 || < 1, then SADM can Бе unconditionally convergent. 


This corollary indicates the reason for which the standard ADM fails in many practical applications (an analogous comment 
is also applicable to HAM Ph. This justifies the need for constructing SADM. Similar to the standard ADM, the solutions by 
SADM are perturbed-analytic; i.e., the solutions are exact near the perturbation origin and a bifurcation appears by moving away 
from this origin. Therefore, SADM is indeed the practical expediency to analytically tackle the nonlinearities. See the example 
we present in Section[3.I]to illustrate this feature. 

In this paper, two symbolic сойе меге developed to conveniently calculate the Adomian terms. 


3.1 | Elucidating Example 


The laminar jet flow over a stationary wall obeys the celebrated Glauert equation?l This third order nonlinear differential 
equation can be formulated in a form resembling the system dynamics, 


d Wa x2 

dt 

d =x) 

at (28) 
=. х) = 0,9) _ 7422 


xO(0) = x@0)=0, x0) = =. ава 


In above, x, x? and х0) are the state variables. This nonlinear equation has а fast exponentially decaying solution at large 
t and is highly resistive to admit an analytic series solution. In this respect, the standard ADM is strictly divergent for Eq. 
(28). Furthermore, the solution provided by HAM, after a careful choice for the auxiliary linear operator, shows a very small 
admissible radius of convergence in the associated h-curve8l However, SADM reveals itself as a successful remedy. To give a 
demonstration of SADM, the second order analytic solution is provided as, 


a) со 2b-aca 
t) =a+bt+ = – ——t 
хо (0 " 2 6 
Q) 2b? +ас » , 2ab? -- aà?c — 5bc 4 
f) =at cet- ———P р 
xo 0) a+c 5 6 
Du Ж. 3 2р2 _ _ 1253 2 
x(t) = c — (25° + ас): + 25 +o Sbe _ ас + 2a°b пае 12Ь? + 5c p (29) 
1 2) 3 
а= xa (tji). b- xa (t1): c= xa (tji). j=1,..,N 
(1),0),63) — 9.2.6) e 
Xp, (15:3) = хы" (5:3), JS2,N 


For the sake of clarity, we define an error for the jet component, x, as, 


N 
„(> / 0a) = + f sea 


t 
И о, 0 


R= (30) 


Roughly speaking, HAM is distinguished by allowing other types of linear operators, or equivalently Es appear. In this respect, HAM may be interpreted in such 
a way that it searches for a contraction through manipulating £o 

3Symbolic code | is associated with a global analytic solution in terms of recursive Jacobian compositions of the state variables for autonomous system dynamics 
with a piece-wise constant control vector (See Section[4]. Symbolic code 2 (which is a standard form for general non-autonomous systems) employs Eq. (13) in a fast track 


mode, that is to obtain А, based on x(t) = x9 + px, +... + p"x,, due to the fact that x,,, does not appear in .A,. 
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In above, 
x9) = oe (t)+ a (0) + ж, (+... (31) 
Moreover, X(t) is the closed form solution for the jet) Table|1|shows a comparison between the analytic series solutions of 
varying perturbed fields with т, = 20 which is an excellent upper bound to fully capture the behavior of the outer flow. The 
computational time for the symbolic codes for different orders of the perturbed fields is tabulated in Table] 
It is worth mentioning that in the zeroth order SADM (P = 0), there is no consideration of the perturbed fields and it is 
straightforward to show that this scenario is equivalent to the explicit Euler method at discrete points. 


TABLE 1 the quality of the residual error R, Eq. (80), for ће elucidating example 


N P P P P P P P 
0 1 2 3 4 5 6 

15 55x10% 17x10? 14x10? 89x10 23x10? 64x10^ 1.3x 1074 

20 538x107 11x10? 18x10? 11x10? 23.5x10^^ 79x10? 7.6x10-9 

30 12x1095 25x10? 8.5х10-* 25x10^ 67x10? 23x10? 2.1x10~° 

40 47x10! $81x10-^ 46x10 58x10 23x10 296x10-5 89x 1077 


TABLE 2 the symbolic computational time for the elucidating example in different orders of the perturbed fields, P 


P symbolic code2 symbolic code 1 


(seconds) (seconds) 
0 0.16 0.03 
1 0.50 0.04 
2 0.68 0.05 
3 0.81 0.05 
4 0.97 0.06 
5 1.15 0.07 
6 1.23 0.09 


Therefore, compared to the existing analytic nonlinear solvers, SADM seems satisfactory to analytically solve the system 
dynamics with a good level of accuracy. 


4 | PERTURBED-ANALYTIC DIRECT TRANSCRIPTION FOR OPTIMAL CONTROL 
(PADOC) 


Without loss of generality, we state the OCP in Mayer form. For simplification, the state constraints are ignored, although the 
direct optimization can easily handle them too. Moreover, the control vector is assumed to be piecewise constant which facilitates 
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understanding of PADOC. 
min J = D(x) tr) 
5.1. 
Tx zx = f (20,00) 
güe) <0 
Po (X(t), to) = 0 
5, (50), tr) 20 

Where, X € RF is the state vector, й € К“ is the control vector, фу € R^, P 


The first challenge is to solve the system dynamics. SADM assumes segments/intervals such that, О = U’ 


following notations should be clarified beforehand , 


2 1 2 L T 
Xo, = (xg O. x D. — xp O) 
C d J 
E 1 2 L 
Xi, = (x19, (0, хо, Os хе O) 
e cT 
*ig, = Хо 


- ($r a) „Èr i Ma Ee T 


> _ (0) (2) (K) T 
Шо, = (ug (D. Ug (D, on (t)) 


Ao (D 4Q Q)T 
Aj = (AO, AD. aO) 


J (885) = 


T. ; 
In above, ( ° ) is the matrix transpose. 
The solution for each element of the state vector is defined as, 


х0 o, = = 2 + xia, FFX 


In above, the analytic-perturbed solutions are obtained ы 


Following the decomposition method, we now proceed with computing Adomian terms of any order. 


< Land 9, € во, О < 
conditions and g € КО is the vector of inequality constraints with 2 having full rank. 


(32) 


L are boundary 


Q, (t, ,. tj). The 


33) 


(34) 


35) 


(36) 


(37) 


(38) 


(39) 


40) 


(41) 


(42) 
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Adomian term of the zeroth order is computed as, 


Ay = f P (ойо) (43) 
Therefore, the Adomian vector of the zeroth order is, 
A, = (адо) (44) 


For the Adomian term of (ће first order, we have, 


0. d () ( LO) (1) x (2) () (0 en 
A = dp А (х Хо, Рта е орд Чанор Риа +. sig J| cs 


@ (45) 
1 
>A = [F(x Xo so, Jas 
Therefore, 
of? af af 
(D — а) (2) (L) 
Ai = E Tm Хо, + а 9) Хо, +... + 3 т) хо, (46) 
хоо, хоо, хоо, 
Eventually, we can write the Adomain vector of the first order as, 
т of . 
A, ==—*19, (47) 
дХо,о, 
Proceeding, as before, for the Adomian vector of the second order, 
af (x? ; üg ) 
> _ 1 4? pe. _14 о? 9,/ Д (ру 
dcm Tap V (Ха бо,)]„— = 21251 a(x? ü ap 58) Ls (48) 
о, 
In an expanded form, it is equivalent to, 
7( zP P(zPo.n B 
Я Е È of (x jio) d? " ү 1 y! д OF (XG ii.) d У рх" 
27 19 ZA 2V0, ! o у. (0) =p NP iQ; 
^ ep) ae MED дух (s) PS (49) 
d (Lp \T 
ap (*0,) |, 0 
After simplification, we obtain. 
> L 
Fi of . 1 д xp 2 
Ay == tn, ts, | o lži, (50) 
Xo, ; 21 2 axon, д 
Following an analogous procedure, for the third order Adomian vector, we get, 
д] 
А,=2х2! x? Dog +3! x 
3 LE UR o]* Жоо, ox n 3.9 
Е (51) 
D D 3 If х9 
+ у ee ax”, с " xio Yio, + ЕИ Ioa, o lži, о, 


Theorem 3. for the system dynamics fx - f (300 10) with the initial boundary conditions Po (Žlto), to) = 0, and u(t) = 


ШУ lig, = ШО, (ues up = 07; (= 5 e SADM over Q, : Q = UŠ Q, (1, ,.t;) is reduced to the following 
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recursive Jacobian compositions' form , 


_ _ пэ _ _ (t = 7*1 dl 
ї) = ї._ + Р, un , М БЕУ GN EE 
Xo) = X(t,-1) > (01) о) (2+1)! 
(1 E us 
= Xt, 4) + Fox, a ) (t7 t1) € + Fy (00), б) ЯЛТ 


X : 4 T : А 
F, is the global perturbed vector of i^" order and со = (с PC ) is the piecewise constant control vector. 
7 j j 


The proof to Theorem 3 is included in an appendix. 


5 | PRIMAL-DUAL COMPLETENESS OF PADOC 


(52) 


Here, we prove that the multipliers associated with PADOC in NLP correspond to a backward analytic solution of the continuous 
adjoint/costate equations in PMP. In addition, it is substantiated that the associated cost functional, in Bolza form, is transcribed 


into an exact analytic form. 
For the OCP defined in[32] the augmented cost is, 
Vi 
7 2 T> | -T> ? IET аз тед 
J = Ф(Х(@„),1,) + V5 Go + Фу «f ато ( FEO. a(t) – 55x) + i | dt 


to 


(53) 


In above, A(t), costate functions, and ji(t) are the time-dependent multipliers, v; and Vy are constant multipliers to involve the 


initial and final boundary constraints. 
The augmented Hamiltonian associated with Eq. is, 


H (X, A) = A f + i^g 
The necessary optimality conditions are, 
OH _ 
ди 
The costate dynamics together with the transversality conditions are, 


0, HT E = 0, azo 


dir | 0H 

dt x 
> дфо T. 
ito) = -| 201] % 

0 
до 28 LIS 
Г REA LEA) 

_ [99 | ;r99; 
Ишден + 59] 


(54) 


(55) 


(56) 


(57) 


(58) 


(59) 
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The cost functional is discretized at the nodes where the analytic solution is partitioned. Therefore, we write the discretized 
version of the cost functional (the Lagrangian) as, 


N 
F = (Xty) ty) + Bo + Фф, + У АТА C (ut) (60) 
j=l 
where A j are the nodal state constraints, 
Age Xo (0) = 0) (61) 
Hence, we have, 
co (1 -t у" 
eg Е эз > MT te 
Aj = (11) - x(t) + 2, (30-0: “Gr (62) 
By defining, л = t; — t; , the above expression is equivalent to, 
027 NX RA(g Qf. At 
Ry = iait У А060) сту (63) 


In above, F, (X, йо) is computed by theorem (3). 

The KKT first-order necessary conditions for the parametric OCP are, 
OX, Qu, Oty (64) 
ÄT =0, >0, k=1,..,.N 


For the interior points, i.e., k = 1,..., N — 1 we have, 


N 
SJ ey ur eg (65) 


k OX А дх = 0 (66) 
From Eq. (63). we can obtain, 

oA 
—*=-] (67) 

OX, 

дА = доз р"! 

——— = EPI [x , 68 
OX, I+ ie йкы) Ga Dl (68) 


In above, I is the identity matrix. Substituting Eq. (67) and Eq. (68) into Eq. (66). 


ð = h”+! 
= ХЕ. + АТ — F, (3p üg) — = 0 69 
k+1 BÈ OX, ( k e) a Ty ( ) 
Next, we show that Eq. corresponds to a backward analytic solution of the system of the adjoint equations by SADM. 
The system of adjoint equations can be expressed as: 
И ðf (X(t), Et 
dap. то А00) (70) 
dt Ox 
Similar to the state equations, the analytic solution for the adjoint equations is developed over an arbitrary interval Q,. By 
choosing £5! = d (+ )dt, we get , 
T of (xq), e(t)) | (71) 


oT _ 3T -1|7 
Дт (0) = а) + £5 È = 
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By expansion, 
7T — 3T 7T 7T 
Ao, (0 = оо, + Ао, Tec Ао, (72) 
Therefore, we can write, 
зт Т 
hoo, = 4 (@,) 


7T _ p-l ӘТ 
Ао, mi Lg, Ву 


(73) 
Т -1 BT 
Ало, = £o, B, 
Where, Br is the transpose Adomian vector of the n” order for the system of adjoint equations. The zeroth term is, 
йы aa) (74) 
0 0,2; ax 
From above, we can note, 
Bj = Mo Ja scs, = 4 (4, ss, (75) 
Therefore, 
Ао = £g Bj = A (JI Css (tj - f) (76) 
For ВТ we write, 
д f (XP ü ) 
ar 4 [үт ?T n7T 9, 9; 
ЙТ = ds | Ma, + pila, sc Pra) |, (77) 
9; 
By simplification, we obtain, 
1; = = 
dr ð Of q т 9f 
ВТ = Д | oe ]«a = (78) 
! 0,9; 2 da Xoo, 19, LO; хоо, 
Р net x x n Ё (D TT n 
We note that the integration is backward; therefore, by substituting хо, апа Ао, їп Ед. (тв) we get, 
L = = = 
Ee д of ү of of 
Br = Bal 2.8 on: + — —— J(u 7! (79) 
| A 2 OX дхоо, (коой) дхоо, дХо,о, Mt) 
Thus, 
Br =A) ITA scs (tj t) = A) IF scs, (t, 1) (80) 
Having Br found, 
2 
t;,—t 
`T _p-lp _ ( j ) `T т, 81 
Ag, ( = La Bi = 21 А ЈАР), (81) 
By induction, we obtain, 
7T 3 (t; = r)” 7T 2 (82) 
до (0 = Lo В,_1 = n! A (t Ja U^, ses, 
Therefore, we can set down the Global Analytic solution for the system of adjoint equations as, 
- 5 а noo "HN (t; = ne m 
ia O = AM) + A, 2, Ј (Р, ба) Ета (83) 
Therefore, in the nodal constraints, we can write, 
>, E =, Е = р"! 
T | jT T 22 = 
А + Agua t А > (Р(Х, ба) ТТ 0 (84) 


п=0 
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Therefore, we can conclude that the multipliers (associated with nodal state constraints) in the direct optimization by PADOC 
correspond to a backward analytic solution of the adjoint equations. The interesting point is that the discrete multipliers by the 
direct optimization can be propagated analytically through the entire interval Q, by Eq. (83). For the boundary conditions, we 
proceed, as before, by differentiation with respect to X(t y), 


OO(X(ty),ty) rO (Rtn) л) 1. > 

a + | || -i0 (85) 
OX(t у) OX(t у) 

This is equivalent to the boundary conditions for the adjoint system. 

For the optimality condition, we have , 


ag TT дА, 2T 0g (8) 


= 4 + 
Qu, “ou, * di, (86) 
jeg, =0, A0 k-L.,.N 
Therefore, 
2 Og(u,) = д 2 р"! 
=T k T 2 > 
+A — F (X, d) ———— = 0 87 
Pk ag, t Zoos, idi) e ny on 


The next step is to show that Eq. corresponds to the optimality condition defined by the discretized Hamiltonian, 


OH (Xp, Адый) _ or of гай 


= + = 0, Те 0, iu > 0 (88) 
дй, k дй, Hk дй, Hi $k = Hk 
To this end, we write, 
af ð d- wv 5 (t-ta) 
= x(t; ma — (89) 
dui, oü dt ou, 4 ale (3-9. £a) Fy 
Dividing Eq. by А, 
zT эйе 
И, д2(и 2 2 
PEREUE ару (90) 
h Qu, Qu, 
Where, 
Ik 
- 1 LIAE 
f= LS Pe) (91) 
1 
Therefore, f is replaced with a pure analytic average over the interval Q, which is the most accurate scenario. Moreover, 
ZT 
Hy >T 2 
^ eg (92) 


The equivalence of the transversality condition follows an analogous fashion. 


: : A pil) ; 2 pO) = po = pe) 
Theorem 4. assume that for a specific OCP, тіп. J? = min; j^ . Therefore, A" = AP”. 
Pp ly ü k k 


2j _ Bf _ 3$ _0 
OX, * дщ ! | Oty 

POTE =0, 20 (93) 
k=1,...,N 
Т8, = 0, fm >20 

PO +4 =1,..,N (94) 
т = 1,...‚ М 
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" Я , x 4 2. ра) d 
Proof: Under the assumption of the uniqueness of the solution for the specific OCP, we must have Ar “aJe Eg 


ЕЯ 
= = 2) = = = k k 
Я A pil) . 2 pO) : pO р z Зыр ZpO 
min; J P" = min; J^ the following holds, АР = A. = Ar. Therefore, АР = АР ж 


Corollary 2. assume that theorem (4) holds for a specific unconstrained OCP with M = 1. Therefore, from Eq. (90), 


T2 p. 
2 fF =0 
ir) pL 
p:J^x/-9 (95) 
AT 9 p. 
it 2 f =0 
2 В, В, z, à zl 
РО: (о) Ff =0 (96) 
From Eq. (95) and Eq. (96) 
2 2 Ü m eg <=. Sp Oo s 3. д = 
M +А нА) АВА. А 
( 1 2 wat an^ 235 Noa? (97) 


This means that — 2 fd is the average of j = T with the weights АГ, Therefore, it is reasonable to look at it as an average 
NonLinear Реа (амі Р). 


СогоПагу 3. we must have М = M beforehand, to instate theorem (4). However, in practice, as M increases, min; A 7 po 
approaches drastically to min; „2 P? ‘Therefore, for а specific OCP, there exists M such that theorem ( 4] holds with an extreme 
accuracy. 


Remark 1. for regular OCPs, aNLP leads to a quite sparse nonlinear optimization. More specifically, with М << N, PADOC 
gives multipliers that can be propagated in a series form over the sub-intervals. Next, by positioning multipliers over the nodal 
state constraints (whose multipliers appear in the optimality condition, a = 0), the staircase optimal solutions from aNLP can 
be converted to the distributed optimal solutions. 


In several examples, we show that even M = 1 provides a very good estimation of the continuous control law. 


6 | CONSTRUCTION OF THE PADOC ALGORITHM 


In this paper, the key elements of PADOC are SADM and aNLP. In this section, we present the PADOC algorithm, launched 
upon SADM and aNLP. 

In the PADOC algorithm, the optimal structure is termed regular if the optimal control does not contain singular/bang 
extremals. Otherwise, the optimal structure is termed irregular. 

The PADOC algorithm is designed to handle both regular and irregular OCPs. In short, first, the general optimal structure in 
the form of staircase function is detected. Next, the outcome is fed to either PADOC handling regular OCP or PADOC handling 
irregular OCP. 

PADOC handling regular OCP exploits aNLP theorem апа remark [1). 

PADOC handling irregular OCP exploits a discontinuity function to identify the number of discontinuities as additional 
decision variables. In the case that Hamiltonian is linear on the control vector, the singular extremal is expressed by a series of 
a truncated order in PADOC handling irregular OCP. 
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Now, we present PADOC algorithm for a generic OCP in Bolza form. A generic Bolza OCP can be pictured asl 
t7 
min J = Ф(х(@), x(t 5). totr) + | £(x(t), u(t), t)dt 


to 


S.t. (98) 


Tx = f (x@, u(t), t) 


g(u@), x(t), t) <0 
ф(х(0), X(t p), totr) <0 


The scalar f € [fo, tp] is time, x € wie (Ito tr], R”) is the state vector and u € L® (Ito, trl, R^) is the control vector. J is 
the cost encompassing a boundary objective Ф : R?"*? — R and an integral objective £ : R” x R? x [to, f f] > К. Dynamics, 
path/control and boundary constraints are f : R” x R? х [ty,f,] ^ R” and g : Rx R^x[tg,1,] > R* and g : Ңң?" > R” 
respectively. 

Next, we employ PADOC to analytically transcribe the continuous OCP into a NLP. It should be pointed out that our analysis 
in the previous section is qualitatively portable to the present section. Therefore, for the sake of conciseness, presenting the 
parametric OCP for the above Bolza form is omitted. The only remaining point is that the integral functional is computed 
analytically by PADOC. 

Eventually, if Mangasarian Fromovitz Constraint Qualification (MFCQ) holds at the local optimal solution, then, there exist 
unique multipliers that satisfy KKT conditions and the solution is reachable through parameter optimization. 


6.1 | PADOC algorithm 
Step 1: the control vector is assumed as piecewise constant; i.e., u(t) = О, uo =U 


qe) over Q : Q= И Өд тд. 
Step 2: we solve the state equations over Q; : Q = UN, Q;(t; ,,1), N > M with a truncated order of the perturbed fields. 
It should be pointed out that for autonomous systems, the Adomian terms are computed by symbolic code 1 (associated with 
theoremB) and for the non-autonomous systems we apply symbolic code 2. 
Step 3: we solve the resulting NLP and extract the staircase solution. In the case of regular OCP, we go to Step 4; otherwise, to 
Step 9. 
Step 4 (PADOC handling regular OCP): applying aNLP (theorem апа remark[1) we extract the optimal solution over О, by 
turning the staircase solution to a semi-continuous one. 
Step 5: The optimal solution due to Step 4 is applied as an initial guess over О, and the resulting NLP is solved. 

The following steps (Step 6, Step 7 and Step 8) are designed to extract analytic expressions for the state and control variables 
over the minimum possible number of intervals. 
Step 6: we find the Лр adaptation by relaxing N and evaluating the global squared residual error for the system dynamics defined 


M (c, e (2) 
j=l Co, a, 


as, 
R= Y RO 
2 
i 2 
к= XL E P xo, — (хай) | 4t (99) 


1 
t)dt 
f; — ti ү, / 
Q; 


We start with N := 1 and increase the order of the perturbed fields, P , to see if the residual error meets a predefined tolerance, 
i.e., R < e. If not, the sequence is repeated for N = 2,3, ---. The outcome is N* and P*. Therefore, the intervals, for both states 
and controls, are modified as: Q, : Q = Ur Qu, (t, 4. fy). 


= 
I 
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Step 7: we express the control vector in a series of a truncated order defined as, 


Ug (f) = Y Y Cale (100) 


n=0 m=0 
The above series is the resulting base functions of an exponential-type linear operator in HAM with a rapid rate of convergence 
(ѕее26] and?) 
Step 8: the resulting NLP, initialized by the solution from Step 5, is solved and the analytic expressions for the state and control 
vectors are extracted. 
Step 9 (PADOC handling irregular OCP): we identify the possible discontinuities by the jump function [u](t) = u(t*) — u(t^) 
(for a better resolution of the jump function, we can employ the methods in2l and3}. 
Step 10: it is important to note beforehand that the SADM solution for the system dynamics ensures that the state vector does 
not violate the possible state constraints in any time instants within the entire time domain. The same argument should hold for 
the control vector. To this end, for an OCP with constraints on the control we can either, 


1. Express the control vector in terms of Bernstein polynomials which satisfy the property of a convex hull, i.e., conu(x,,) := 
{ Y aded 4 x0 A LAS 1). Sec?! for more information. 


о=0 0 
2. Treat the number of possible discontinuities as extra decision variables in the frame of NLP and express the control, 
outside the bang region, by Eq. (100). 


Here, we have followed the latter scenario. For K number of discontinuities, the additional decision variables are f. k-l.,K 
and the time domain is partitioned as, 


K 
9 - [ Jo, (7, 7) (101) 


К=1 


Step 11: hp adaptation is checked over the subdomain Q,. For each extremal Q,, SADM is applied over Af,, 


e k-1 k 

ye 1 

Ai, = "Ww ae de 1+ Nau 9 Na BS Tp K (102) 
i=l i=l 


The global residual error can be rewritten as, 


же c fi risa am 


We start with N, :— 1 and increase the order of perturbed fields and check if R, < e. If not, we set N, = 2,3, ... and repeat 
this step. The outcome is Рё and NT. 

Step 12: over the bang extremals (Q5), the control vector is set as piecewise constant. As for the singular extremals (Оу), we 
use the series function, Eq. (100). up to a truncated order. 

Step 13: we solve the resulting NLP initialized by the staircase solution in Step 3. 


7 | ILLUSTRATIVE EXAMPLES 


This section is devoted to some clarifying examples. For the sake of brevity, PADOC is inspected in more details throughout 


the first four examples [7.3] [.4). However, for the next examples [7.6), we have mainly displayed the quality of 
the optimal solutions by PADOC and comparisons wherever applicable. 
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In this research, interior-point/barrier algorithm of the NLP solver fmincon from MATLAB® Optimization Toolbox was 
adopted as the optimization module. The codes within the present paper were developed in MATLAB®, R2020b, on a laptop, 
configured by Intel(R) Core(TM) i7 (10TH GEN)-10750H CPU @ 2.60GHz. 


7.1 | Orbital Transfer 


The OCP of transferring a rocket from an initial orbit to a target orbit in a fixed time is expressed а 
: 1,2 2 1 

= Ф(ї,) = – | – t) + t — 

min = (1) = – |5 (^0 ) + (ty) 129) 


S.t. 


U 
d y (104) 
—u- дар) 
r2 


os au + Acos e(t) 
r 


А = 0.01, (0) = 1.1, 0(0-20, wu0)20, vO)= =, t; =50 


V1.1 
The control variable is the thrust steering angle measured from the local horizontal e(t). 

Interestingly, we can show, by order of magnitude analysis, that the solution to Eq. is in close proximity to the solution 
with e(t) = 0. Therefore, the objective is almost unaffected by the control (to be compared later). Moreover, in a mathematical 
point of view, this is a triple state OCP since 0(7) is isolated. 

For this OCP, Hamiltonian can be written as, 


2 
M Au A( -5 + Asine) +4,(- “© Асос) (105) 
F r r 


Therefore, PMP for the optimality condition, Hi = 0, gives, 


1 40) 


e(t) = tan^ (106) 
AQ) 
SADM over О j reads, 
= Asi 1 v? 2 
PO) = ria, +U, tt (Asineg — ae dun m 
ps Ad 1 1? ; 
u(t) — Wer t ( sin €o, — E T a 
2 
€ (5 - 7 e P (Asinen, - “©)) P+... (107) 
rà r2 r Р 
: ир 
v(t) = 04 + (А sin Eg, — un 
+ (2 - *(Asine EN үт -*)) P+ 
Dor G p'r r 9; rm u 


1-1 

For this OCP, aNLP can be initiated by setting М = 1. This is equivalent to e(t) = ё. SADM of (ће 1*' order was confirmed 
to be quite accurate and stable over 50 nodes. Therefore, we take N — 50. 

From Eq. (106), the optimal control is a function of 4, (1) and 4, (7). Therefore, it is reasonable to ignore positioning the nodal 
state constraints for A,(t) in aNLP. 

For 50 nodes, 49 segments, and M = 1, the size of the aNLP is 99 variables. Next, we solve the aNLP by the optimization 
module and extract multipliers À,, and 4,,. Eventually, we plug these multipliers into Eq. and the optimal control law with 
the resolution of the state nodes (50 nodes) is extracted. Fig. Ша shows the performance of aNLP with M = 1,5,10. For the 
initial guess, we used a constant steering of 0.001 rad. 
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Clearly, aNLP, even with M = 1, provides a very good topological approximation of the continuous control. Table (3) displays 
the performance of the optimization module on independent aNLP solutions. 


TABLE 3 Performance of the optimization module on aNLP in example 


M iterations Optimization Time Optimal Cost 


(seconds) 
1 3 0.17 0.095841 
2 7 0.29 0.095734 
5 12 0.38 0.095184 
10 16 0.72 0.095132 
50 44 3.2 0.095118 


In fact, for M = 1, the actual aNLP has only one decision variable and the existence of the nodal state constraints is not 
essential and NLP can indirectly retain the concept of the multipliers associated with the nodal state constraints. In this respect, 
for the aNLP with M = 1, we have performed a direct optimization over the only existent decision variable, e(t) = e. The main 
aspect of this case is to provide a qualitative visual on the influence of different transcription methods on the Hessian search 
direction, accuracy and the effective domain for the initial guesses. Fig. [i}b shows the behavior of the objective as a function 
of e(t) = Є in applying second order SADM in the PADOC algorithm. Fig. Ш shows an analogous behavior in а stable Euler 
method. 

In the end, in the case that theorem (4) holds, there seems to be a complementary theorem that relates the topological local 
optima from the univariate ring (Fig.[I}b and[I]c) to a multivariate function. 

The optimal states from the PADOC algorithm are shown in Fig.[I}d and[I]e where we can verify that the optimal solution is 
in close proximity to the solution without having a control. For the case of e(t) = 0, we find Ф(7 p — 0.096247. 

Table[H]shows the global squared residual error (Eq.[99) for SADM of a varying order to manifest the hp-adaptive quality for 
this OCP. 


TABLE 4 the A p-adaptive quality, Eq. (99). for example 


N P P P P P P 

0 1 2 3 4 5 
20 38x10 3.6х10-! 297x10^^ 24х10% 84x107 9296x107 
30 2.5x10P 58x10? 29x10* 65x10° 12x10? 33x10 
50 38x10? 86x10 76x10? 9x10 3x7 2x 1078 
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7.2 | Van der Pol Oscillator 


The problem is defined а$541 
2 


min J = Ju (EAO + x0) + u(t) dt 
0 
S.t 
d (108) 
a = х 
La = —x; + x,(l — x) + u(t) 
x,(0)=1, x,(0)=0 
The Hamiltonian for this problem is, 
H = TG + E t uw) + A, х2 + Ax, ( —xX,+x(1- х) + и) (109) 
Therefore, form РМР (= = 0), 
u(t) = —A, (t) (110) 


For this OCP, we used SADM of the first order over 40 nodal state constrains (N = 40) with M = 1,2,5. The multipliers 
were extracted from the optimization module and together with Eq. (110), the optimal control was found (see Fig. [2}a). Indeed, 
a very good agreement is verified. It should be mentioned that the nodal state constraints were distributed over x,(t) since the 
optimal control has no dependency on 4, . Table|5|shows the hp-adaptive quality for SADM (Eq. 9). Following the PADOC 
algorithm (Step 6-8), we were able to extract high order optimal series solutions for the state equations as well as highly accurate 
piecewise interpolants for the optimal control. The minimum cost by PADOC is J = 1.04285. Table [6] shows a comparison 
between the semi-analytic solution reported in®4l and PADOC. Рів. exhibits ће optimal series solutions by PADOC for the 
state functions. Table[7|shows the performance of the optimization module on aNLP with M = 1, 2, 5. 


TABLE 5 the Ap-adaptive quality, Eq. (99), for example|7.2 


N P P P P P 
0 1 2 3 4 

2 37х10! 32x10? 16x10? 26x10^ 1.2x107 

3 35x10! L1x10? 65x10^ 13x10^* 34x10? 

4 31x10! 91x10^ 3.1x10 54x10? 86x10 

5 23x10! 63x10^ 18x10* 25x10^?- 5x107 


TABLE 6 comparison between PADOC and Pade-HPM™! for the optimal control u(t) in example|7.2 


t PADOC  Pade-HPM 
0 0.010914 -0.011719 
0.4 0.517903 0.514115 
0.8 0.721478 0.720252 

1.2 0.644546 0.645724 

1.6 0.370995 0.372690 

2 0.009343 0.000775 
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TABLE 7 Performance of the optimization module on aNLP in example 


M iterations Optimization Time Optimal Cost 


(seconds) 
1 4 0.12 1.10339 
2 8 0.19 1.10205 
5] 18 0.32 1.05747 


7.3 | Brachistochrone Problem 
This is among the first appearing OCPs and was solved by Johann Bernoulli by calculus of variations. Brachistochrone OCP 
can be defined as3] 


min J = t 
u(t) J f 


5.1. 
а р 
qi = vsinu(t) 


: (111) 
Z Ag me t 
d cos u(t) 


d 
—p= t 
230 g cos u(t) 


x(0) = x(0) = v0) 20, x@;)=yt;)=2 
Let us first consider M = 1. The solution of the system dynamics by SADM gives the exact closed form solution. In other 
words, since the higher order perturbed fields become zero in the Adomian recursion scheme, the induction of the close form 
solution is immediate. Therefore, the closed form solution by SADM reads , 


x(t) = 2 sino” 
y(t) = 5 cos ^r" (112) 
v(t) = g cos(a)t 


For M = 1, the corresponding aNLP, can be put into the following form, 


s а PO) m ~ (8 Е 
Д=ї,+%, G sin(2a)t, — 2) + 05 cos^(u)r^. — 2) (113) 
Therefore, the optimization module attempts to solve, 

o =1+0; : sinQi)t , + 0, g cos (ü)r , = 0 
0j . 8 LAU. ud uuu us 
Ou = 53 cosQü)t", — us sinQüjr, =0 (114) 
$ 5 —42 = 
4 sinut, -2 =0 
& 


20-32 = 
z 008 (шг, —2=0 
The close form solution reads , 
8 . 0. 9 2 


л 
П УЕНЫ ee ES (115) 


и = 


The Ap adaptation was found with R = 0 upon Р = 1 and N = 1. In addition, it was found that u(t) = at best describes the 
exact optimal control. 
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On processing the SADM, we found that Adomian terms of the second order forward are all zero. Therefore, the SADM of 
the 1% order is the exact close form solution. Consequently, given u(t) = at, the close form solution by SADM becomes, 


t  sin(2at) 
Me G i i ) 
= Č sin? 11 
y(t) = 525 sin'(ar) (116) 


v(t) = © sin(at) 
a 
In this sense, the augmented cost is, 


8f,  gsinQat,) А g ту 
L eg —2) + (E sin an -2) (117) 
In this situation, the solution by the optimization module for t , and a corresponds to the solution of the following trigonometric 
algebraic equations, 
8t,  gsinQat,) 2 
2a 4 (118) 
$ u2 = 
23 sin'(at;) - 2 = 0 


The close form solution reads, 
(119) 


In numerics, the optimal solutions are a — 1.4629977 and t, — 0.8243387 which is comparable to the aNLP result (with 
M = 1), t; = 0.9030473. 


7.4 | Aly Problem 


This linear-quadratic regulator OCP is expressed a2] 
5 


. 1 2 2 
mg = ; Гао + 20да 


0 
5.1. 
(120) 
qu => Xa 
d 
T = u(t) 
|u(t)| x 1 


Through the PADOC algorithm, this problem was solved over 2 segments; (0, t4) and (t4, 5). t4 is the identified discontinuity. 
We found that the Adomian terms of the second order forward were all zero. Therefore, the optimal solution by PADOC for this 
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OCP is in close form. The optimal solution by PADOC is as follows , 
t € (0,14) 
u(t) = —1 
x(t) = siga 
Pes 
x,(t) 21—t 
t € (tz, 5) 
u(t) = ae™ +b 


(121) 


х|@ = ae + 22 + (1-1, — bta ae) t4 — xf 


х›(0) = —ае + bt -- ae — bt, —ty +1 


The decision variables were optimally identified as; t; & 1.405714, a ғ 1.495751, b ғ 0.012781. Eventually, the optimal 
cost was computed as, J ~ 0.364513. The switching time is comparable to that inl ty #5 1.413764. The optimal solutions are 


graphed in Fig. Ba and[3}b. 


7.5 | Maximum Range of a Hang Glider 
This OCP is stated аѕ29} 


min J = —x(t 
m (ty) 


—x =D 
dt С 

Г 

d ° 

а 1 . 

dit - -( — Lsin(y) — D cos(n) ) 

Lo, = -(u cos(y) — D sin(y) — w) 

O<c,(t) X 1.4 

х(0) = 0, у(0) = 1000, 0,00 = 13.23, v, (0) = —1.288 (122) 


уі) = 900, v(t) = 13.23, ьт) = —1.288 


п = tan! [Em v, = V v2 + (v, — uo) 


Uy 


L= TA D= 5pcp (ex) Su. W = тв 


х 2 2 
cp(cz) = ср + кс], u(x) = us C06 125) (1 = is = 25)?) 
и =2.5, К= 100, сь = 0.034, К = 0.069662, т = 100 


а,тах 


S=14, р= 1.13, в = 9.81 


The values аге given іп SI unit. The authors inP9 have reported that the convergence cannot be reached for the above OCP in 
its complete form upon applying a direct collocation method of a cubic Lobatto type. For this reason, they combined direct and 
indirect approaches together with parameter continuation technique to solve this OCP. On trying PADOC, The Ap adaptation 
was reached for 1° order (above about 120 nodes), 2"4 order ( above about 80 nodes) and 3'7 order (above about 50 nodes). It 
is noteworthy mentioning that there was no sensitivity by PADOC to an initial guess within a practical range. The maximum 
range computed is, x(f,) = 1247.947903 with the flight duration, f, = 98.409791. The computed quantities agree very well 
with the report in?) that is x(t) = 1247.60 and т, = 98.380. Fig. [а shows a comparison between the optimal control 
by PADOC complete solution compared with the optimal control by aNLP (without interference of the Hamiltonian) with 30 
segments, 150 nodal state constraints and SADM of the second order. Fig. /4}b, с. а апа Ве manifest the optimal states by 
PADOC complete solution. Table[8]shows the time performance of the optimization module with M = 2,3, 10, 30 and second 
order SADM. It is noteworthy mentioning that in principle, the number of the decision variables needs to exceed the number of 
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boundary equations. In the case of equivalence, the aNLP turns into a shooting-like scenario, that is to find the decision variables 
such that the boundary equations are satisfied. For M = 2, the aNLP contains 2 segments for c, (t) plus f +. Since the OCP has 3 
boundary equations, this case coincides with a shooting scenario with no actual decision variables in the sense of optimization. 
The solutions tabulated in Table [8]are independent with a fixed initial guess. We continued the PADOC algorithm taking into 
account the aNLP solution with M — 30 as the initial guess and the final solution was reached in about 26.35. 


TABLE 8 Performance of the optimization module on aNLP in example 


M iterations Optimization Time Flight Duration Optimal Cost 


(seconds) (seconds) (meters) 
5 0.13 90.53366 1199.22208 
3 9 0.26 93.58561 1208.61251 
10 34 0.95 94.80555 1219.70867 
30 69 8.4 98.03995 1244.43072 


7.6 | Minimum Time Aircraft Trajectory in Climbing Phase 


The original system dynamics governing an aircraft in its climbing phase contains variables with different time scales, mass 
m(t) and air slop y(t). A normal remedy is to apply singular perturbation (for more details, sec). Under a singular perturbation 
analysis, the system dynamics becomes as follows. Thrust and fuel flow are modeled by Base of Aircraft DAta (BADA) model. 
In addition, air density is approximated by International Standard Atmospheric (ISA) model. 


qupd = ty 

Lh = vy(t) 

To x ш = Laws (сь, + Cn, (28 .)) ey) 
£m = -C, T0) 

Ymin € Y(t) < 0.262 


h(0)- 3480, v(0) = 151.67, m(0) = 69000 (123) 
h(t) = 9144, (ту) = 191, m(t) = 68100 


һ U 
Gn ee iren С=С (1+ 2 ) 


P(h) = By 


AOV, O(h) = Oy -ph pU - = 
0 


S= 122.6, g=981, Су = 141040, Сү, = 14909.9 
Cp, = 6.997е710 Cp = 0.0242, Cp, = 0.0469, С, = 1.055e? 
C, = 441.54 К = 287.058, Ө, = 288.15, f = 0.0065, Р, = 101325 


This singular OCP has been solved in® for Ymin = —9.262. Here, we have relaxed y,,;, to analyze the performance of PADOC 
in a varying control bound. The PADOC algorithm confirmed a bang-singular-bang solution for this OCP. Next, we partitioned 
the domain by the global residual error as instructed through the PADOC algorithm. The Лр adaptation was confirmed as (with 
the R well below 1075), one segment over (0, 14.) 5 segments over (ta, Я ta) and | segment over (ta, t,), each segment with the 
perturbed fields of the 4^" order. Next, we used the interpolant series, Eq. (100). for the singular arc and solved the resulting 
NLP. For the case of ymin = —0.262, Fig.|5+a shows a comparison between the report іп57 (27 = 644.2), aNLP with М = 50 
(27 = 644.52101) and ће PADOC complete solution (27 = 642.34312). В.Б displays the optimal controls for y,,,,,, = —0.262 
to Ymin = 0.022 with the step, dy,,;, = 0.01. Fig.Ble betrays t% as a function of y,,;, which shows a marginal affectibility. Fig. 
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pja, le and BHE exhibit the optimal solutions for the state variables. Table P]depicts the time performance of the optimization 
module with M = 2,5, 10, 50 and the second order SADM. The final solution was reached in about 8.35 on using the results 
from the aNLP with M = 50 as the initial guess. 


TABLE 9 Performance of the optimization module on aNLP in example 


M iterations Optimization Time Minimum Time 


(seconds) (seconds) 
2 7 0.14 658.10901 
5 28 0.92 651.79756 
10 67 1.8 647.40787 
50 113 6.6 644.52101 


8 | CONCLUSION 


There exist various numerical integration and collocation schemes to solve the system dynamics and to form an NLP and only 
some of these methods are complete. Reviewing the literature, it appears that the existing transcription methods, have been 
introduced as alternatives with limited studies on the aspects of the superiority of one over another. Therefore, the ultimate goal 
will be to introduce PADOC as a unified framework for direct trajectory optimization. 

Nevertheless, there has been no general analytic transcription method in the literature until now. In short, the reason lies behind 
the fact that the area of research belonging to the analytic nonlinear solvers is, in some respects, immature and the so-far known 
analytic methods are successful only on some types of nonlinearity and they do not guarantee convergence for a wide class of 
nonlinear problems. This gap was replenished by introducing a powerful analytic nonlinear solver, SADM. This new analytic 
solver together with aNLP were deployed to form PADOC which is a complete method since it allows an order-preserving map 
between the dualized problem and the discretized problem. 

PADOC was run over many OCPs and its robustness, computational efficiency and reliability were confirmed. Here, some of 
the distinct attributes of PADOC were manifested through the examples. Specifically, Ex[7.1|(an orbital transfer problem) and 
Ех[72] (the Van der Pol oscillator) were mainly devoted to the role of амі Р and how it reduces the decision variables while 
retaining the accuracy. In Ex. [7.3] (Brachistochrone problem) and Ex[7.4] (Aly problem) we mainly showed the potentiality of 
PADOC to extract nearly close form optimal solutions. Eventually, through Ex[7.5](maximum range of a hang glider) and Ex{7.6] 
(minimum time aircraft trajectory in climbing phase) we demonstrated the performance of PADOC in more complex arenas 
with highly accurate outcomes. 

The purpose of this research was mainly to introduce PADOC not only as an algorithm but also as a new prospective method 
to absorb the attention of the research community. 

In view of the present authors, the potential extensions to PADOC will be 1) to scrutinize the sparsity of the search direction 
matrix under PADOC compared to other transcription methods; 2) to explore other linear operators, leading to topologically- 
different perturbed fields and to analyze the CMP criterion; 3) to extend aNLP for constrained and singular OCPs as a general 
bridge between the differential theory and the variational theory through Fixed Point Property (FPP); 4) to compare the optimal 
results by PADOC with other transcription methods in terms of accuracy and the computational time for practical applications. 
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| APPENDIX: PROOF OF THEOREM 3 


Proof: The autonomous system of 1“ order differential equations over an arbitrary interval Q; is integrated with се = J (•)а t 
1 


as, 


j-1 


Xo (0 = Xo 4) + / FEO, ča Jdt (124) 


j-1 
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Following SADM, we assume, 


noo 


5a (0 = Ys, 0,2) = У A. A= lap ол (125) 
On using into Eq. (124). Е 
хоо, = хо (1) 
t 
Xia = J Ай, il (126) 


RI 
St 
м 

£ 

n 


Let us define an auxiliary vector F, (Х, с, 


o(X, Co, ) m £p. бо) 
= 5 X (127) 
f (3.45) = (P aes) Fa) #21 
According to Eq. (126), 
t t 
X19, = / Ада = / f(x. Co) xq, 04 = f (3, Co, ) sx, 0б! т tj) (128) 
tj ti 
= Fo (x). ĉo) (t =t) 
From Eq. (47), we have, 
Ay = Iu PF G.S sas (E ta) = PG E) 7 ti) (129) 
Therefore, 
t t 
о = J А, = / JU (X. ča, ) ssa, C —t)_,)dt 
10-1 lj (130) 
(t-t) 
= Р, А са) pu 
From Eq. and by making use of Eq. (128) and Eq. (130 
> t > 
A, = PEN Gu у us f (S.S) |F(% ta) (131) 
=х@, 1) 
It can be shown that, 
ps + Xs D LPT e а) = JSAP) (132) 
Therefore, 
(1-51) 
Ay = (006) 2 (133) 
апа, 
qe ba) 


X39, = " Аа! = Fa (žlta) бы = е (134) 


ta 
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Following the same procedure for the perturbed vectors of higher orders, we obtain, 


(t-41) 
5 S pate. 
X40, = | Аа = Р, (0,1), čo) — 
tji 
(135) 
t В 
Ши ЧР m3 
Ўр = А Aidt = F (3t), ča) — = 
tja 
Since, Xo, = Xoo, + Xio, + + Хо, 
X Y P(x T BÍ > (rta) 
Xo (t) = X(t;  * Fo(X(; 9.8 )(t— t)-1) + + F (X54. Co, ) ————— M (136) 
і | i" (п+1)! 


Corollary 4. for an autonomous system of 1° order differential equations as defined in theorem, the global recursion scheme 
for Adomain vectors is, 


"um к эш (137) 
n " a (Aa ) Ap 
Proof: From Eq. (135), we can write, 
4 т (> > (t Е tj- ) 
А, = F, (Xo, ё) —— (AR 
Therefore, we have, 
_ " _ _ (t _ fü у! 
А, = Fai (Xoo Co )——— и 
апа, 
PM teta) s a eta 
A, = F, Goa, ёа) се. = Т К а EE p 
From Eq. and Eq. (140). 
2 [—1.: ES =p 
A, = кые лд (141) 


Remark 2. The Global Analytic solution presented in theorem З|сап be extended to a large class of non-autonomous п!” order 
nonlinear ODEs. However, derivation of such a formula is not within the scope of the present work. 
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FIGURE 1 Example Pala: the converted aNLP solutions compared to the PMP solution. (b): the objective in PADOC aNLP 
with M = 1 as a function of control. (c):the objective in Euler-based transcription with M = 1 as a function of control. (d) and 
(e): state variables from e = 0 (red lines), compared to those from PADOC algorithm after optimization (black lines) 
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FIGURE 2 Example[7.2] (a): the converted aNLP solutions compared to PMP solution. (b): the state variables from PADOC 
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FIGURE 3 Example (a): the optimal control from PADOC algorithm. (b): the optimal state variables from PADOC 
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FIGURE 4 Example [7.5] (a): the PADOC aNLP with M = 30 (red line) compared to the PADOC complete solution (black 
line). (b),(c),(d),(e): optimal state variables from the PADOC complete solution 
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FIGURE 5 Example [7.6] (a): PADOC aNLP with M = 50 compared to the PADOC complete solution and PMP combined 
direct method. (b): the PADOC complete solutions in varying у,„;„. (c): minimum time as a function of Yin. (d),(e),(f): the 
optimal state variables in varying ymin 


